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Abstract 

We study Maxwell’s equations in random media with small fluctuations of the electric permittivity. We consider 
a setup where the waves propagate toward a preferred direction, called range. We decompose the electromagnetic 
wave held in transverse electric and transverse magnetic plane waves, called modes, with random amplitudes that 
model cumulative scattering effects in the medium. Their evolution in range is described by a coupled system of 
stochastic differential equations driven by the random fluctuations of the electric permittivity. We analyze the solution 
of this system with the Markov limit theorem and obtain a detailed asymptotic characterization of the electromagnetic 
wave held in the long range limit. In particular, we quantify the loss of coherence of the waves due to scattering by 
calculating the range scales (scattering mean free paths) on which the mean amplitudes of the modes decay. We also 
quantify the energy exchange between the modes, and consequently the loss of polarization induced by scattering, 
by analyzing the Wigner transform (energy density) of the electromagnetic wave held. This analysis involves the 
derivation of transport equations with polarization. We study in detail these equations and connect the results with the 
existing literature in radiative transport and paraxial wave propagation. 
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1. Introduction 

Understanding the interaction of electromagnetic waves with complex media through which they propagate is of 
great importance in applications such as radar imaging and remote sensing [1, 2], optical imaging [3], laser beam 
propagation through the atmosphere [4], and communications [5]. As the microstructure (inhomogeneities) of such 
media cannot be known in detail, we model it as a random held, and thus study Maxwell’s equations in random 
media. The goal is to describe features of the solution, the electromagnetic wave held, which do not depend on 
the particular realization of the random medium, just on its statistics. Of particular interest in applications are the 
statistical expectation of the solution, which describes the coherent part of the waves, and the second moments which 
describe how the waves decorrelate and depolarize, and how energy is transported in the medium. 

In most applications the inhomogeneities are weak scatterers, modeled by small fluctuations of the wave speed. 
They have large cumulative scattering effects at long distances of propagation. Among them are the loss of coher¬ 
ence, manifested mathematically as an exponential decay of the statistical expectation of the wave held, and thus 
enhancement of the random fluctuations, and wave depolarization. The quantification of these effects depends not 
only on the amplitude of the fluctuations of the wave speed, but also on the relation between the basic lengths scales: 
the wavelength, the scale of variations of the medium (correlation length), the spatial support of the source, and the 
distance of propagation. 

Recent mathematical studies of electromagnetic waves in random media are in [6] for layered media, in [7] for 
waveguides, and in [8] for beam propagation in open environments. They decompose the electromagnetic field in 
plane wave components, transverse to a preferred direction of propagation called range, and then analyze the evolution 
of their amplitudes, which are frequency and range dependent random fields. The details and results differ from one 
study to another, because the geometry and the scaling regimes are different. For example, the decomposition in 
waveguides leads to a countable set of waves, whereas in open environments, as we consider here, there is a continuum 
of plane waves. The regime in [8] leads to statistical wave coupling by scattering, however the waves form a narrow 
cone beam and they retain their initial linear polarization. 
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In this paper we build on the results in [7, 8] to obtain a detailed characterization of polarization effects in random 
open environments. We consider a medium with random electric permittivity that fluctuates on a scale (correlation 
length) t that is larger than the wavelength d by a factor l/y, with j 6 (0,1), and the waves propagate over many 
correlation lengths. By assuming that the support of the source is similar to the correlation length and that the 
fluctuations of the medium are small and smooth (i.e. their standard deviation a is small and their power spectral 
density decays fast enough), we identify an interesting regime where the waves propagate along a preferred direction, 
the range. It is between the paraxial regime studied in [8], where the waves travel in the form of a narrow cone beam, 
and the radiative transfer regime in [9, 10], where the waves travel in all directions. In our regime the waves propagate 
in a cone whose opening angle is significantly larger than in the paraxial case, but smaller than 180 degrees, so that 
the backscattered waves can be neglected. The validity of this regime is controlled by the parameter y and the distance 
L of propagation. We take L :s> { - Ajy so that the scattering effects in the medium build up, but to ensure that the 
waves remain forward propagating, we restrict the distance to L ~ 

From the mathematical point of view, the advantage of having a preferred direction of propagation is that we 
can reduce the analysis of Maxwell’s equations to the study of the range evolution of the random amplitudes of the 
components of the wave held. These amplitudes satisfy a system of stochastic differential equations driven by the 
fluctuations of the electric permittivity, and can be analyzed in detail using the Markov limit theorem [11, 12]. Our 
main results are: 

1. The quantification of the scattering mean free paths, the range scales on which the components of the electro¬ 
magnetic wave held lose coherence. 

2. The quantification of the statistical decorrelation of the waves over directions. 

3. The derivation of the transport equations for the energy density, which allows us to quantify the depolarization 
of the waves and the diffusion of wave energy in direction. 

4. The connection of the results to the radiative transport theory with polarization given in [13, 9, 14], and to the 
moment equations associated to the random paraxial wave equation given in [15, 16, 17] (scalar case) and [8] 
(electromagnetic case). 

The paper is organized as follows: We state up front, in Section 2, the mathematical problem and the main results. 
The derivation of these results uses the formulation and scaling described in Section 3. The derivation of the wave 
decomposition is in Section 4. To give an intuitive interpretation of the decomposition we consider first homogeneous 
media. Then we give the decomposition in random media and derive the system of stochastic differential equations 
satisfied by the random wave amplitudes. The Markov limit of the solution of this system is obtained in Section 5. 
We use it in Section 6 to quantify the loss of coherence of the waves. The analysis of the second moments of the 
amplitudes and the derivation of the transport equations is in Section 7. They are connected to the radiative transport 
theory in [13, 9, 10] in Appendix B. We illustrate the results with numerical simulations in statistically isotropic 
media in Section 8. We give a detailed analysis in the high-frequency limit y —> 0 in Section 9 and connect these 
results to the white-noise paraxial wave equation studied in [8] in Appendix C. We end with a summary in Section 
10 . 

2. Statement of results 

We begin in Section 2.1 with Maxwell’s equations satisfied by time-harmonic electric and magnetic fields in 
random media. The scattering regime is defined in Section 2.2, by identifying the important scales and describing 
their relations. The wave decomposition in transverse electric and magnetic modes is stated in Section 2.3. The 
characterization of the first and second moments of the random amplitudes of these modes are our main results, stated 
in Section 2.4. They are proved in Sections 6 and 7. 
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Figure 1: Geometric setup. The source of diameter X is localized at the origin of range z, and emits waves in the direction z, in the random medium 
extending from z = 0 to z = L. 


2.1. Maxwell's equations in random media 

The time-harmonic electric and magnetic fields E{£) and &{£) satisfy Maxwell’s equations 


V X E{x) - iojpoH{x), 

(1) 

V X Hi£) = J{£) - iue{£)E{£), 

(2) 

V-[e(it),E(a?)]=p(^), 

(3) 

V ■ [p„i?(a?)] = 0, 

(4) 


for e and frequency cu e K. The setup is illustrated in Figure 1. The excitation is with the curi'ent source 
density JT localized at z = 0, that emits waves in the direction z, called range. The system of coordinates is a? = {x, z), 
with vector x - (xi,X 2 ) in the cross-range plane. The waves propagate in a linear medium with random electric 
permittivity s{aS) and constant magnetic permeability fig. The medium is isotropic and non-absorbing, meaning that 
e{£) is scalar valued and positive. Extensions to variable magnetic permeability and to dissipative media with complex 
valued permittivity tensors are possible, but for simplicity we do not consider them here. 

We focus attention on equations (1) and (2), because the other two equations are implied by them. Equation (4) 
follows by taking the divergence in (1) and the charge density p{x) in equation (3) is related to the current source 
Jri^) by the continuity of charge equation 

iLopi^) - V ■ Ji^), 


obtained by taking the divergence of equation (2). 

The random electric permittivity models small scale inhomogeneities in the medium, 
constant value as described by 

er(^) = 1 + 1(0,L)(2)«V 



It fluctuates around the 

(5) 


where er(3^) = is the relative electric permittivity, and v is a dimensionless stationary random process of 

dimensionless argument in The process has zero mean 


]E[v(r)] = 0, 


( 6 ) 


and autocorrelation 

E[v(fi)v(f^)] = ■F'), Vr, f^ e K^. (7) 

We assume that v is bounded and differentiable, with bounded derivative almost surely, and that is integrable, 
with Eourier transform (power spectral density) 


JW' 


( 8 ) 
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This is a non-negative function by Bochner’s theorem, and we suppose that it is is either compactly supported in a ball 
in of radius Oi l), or it is negligible outside this ball. The autocorrelation is normalized so that 



= Oil), 


3?(0) = (9(1). 


(9) 


The length scale {in (5) is the correlation length and the positive and small dimensionless parameter a quantifies the 
typical amplitude (standard deviation) of the fluctuations. 

The indicator function 1(o,l)(z) in (5) limits the support of the fluctuations to the range interval z 6 (0, L). This 
allows us to state easily the boundary conditions satisfied by E and A. The truncation at z = L may be understood 
physically in the time domain, where L is determined by the duration of the observation time. For a finite time the 
waves emitted by the source cannot be affected by the medium beyond some range, called here L, which is why we 
can truncate the random fluctuations there without affecting the solution. The truncation at z = 0 is consistent with the 
fact that in our scaling the backscattered waves are negligible. Thus, the waves that propagate from the source in the 
negative range direction have no effect on the waves at z > 0, which is why we can truncate the fluctuations at z = 0. 


2.2. The scattering regime 

There are four important length scales that define the scattering regime; the wavelength A, the correlation length 
the distance of propagation L, and the support X of the source in cross-range. The wavelength is defined at the 
reference wave speed Co = 1/ yjSoiJ^o by T = IjiCoIcj, and L is of the same order as L. The support X of the source in 
cross-range determines the initial opening angle of the cone (beam) of emitted waves. To define it we let the current 
source be of the form 

( 10 ) 

where J(t’) = {J(r), Jz(r)) is a vector-valued function of the dimensionless vector r e The magnitude of S is 
assumed negligible for |r| > 1, so X is the diameter of the spatial support of the source. The scaling by the constant 
impedance is chosen for convenience. 

Our scattering regime is defined by the scale ordering 

A<(~X^l, (11) 

and the small standard deviation a of the random fluctuations of the electric permittivity, satisfying 

~ til. (12) 

Here the symbol ~ denotes the same order. By taking a large distance of propagation L, as in (11), we ensure that the 
cumulative scattering effects in the random medium have a significant effect on the electromagnetic wave field. This 
effect is summarized in Section 2.4 by the exponential decay in z of the statistical expectation of the wave field, which 
is a manifestation of its randomization, the exchange of energy between the wave modes, and the diffusion of energy 
over directions of propagation. Because of this diffusion we restrict the propagation distance L as in (12), so that the 
forward going waves emitted by the source remain forward going for all z e (0, L). 

2.3. Transverse electric and magnetic modes 

The interaction of the waves with the random medium depends on the direction of propagation, which is why 
we decompose the electric and magnetic fields in plane waves with wave vector it. The decomposition stated here is 
derived in Section 4.4, and in our scattering regime it consists of two types of forward going waves which we call 
modes: transverse electric and magnetic. 

The decomposition of the electric field is 

f + a-^(/(,z)il'^(K)]e''''^'^, (13) 

J|q<i ( 2 ^) 


4 





(14) 


and for the magnetic field we have 

H(x) = f ^^/ 3 ^H/c)[a(/c,z)u'^(/c) - a-^(A:,z)H(*:)]e'“'^, 

Jw<i (2^) 

where k = InjA is the wavenumber and we use the notation dikK) - k^dK. The three-dimensional wave vector is 

JI^{k,I3{k)), K^{Ki,K2)eR^. (15) 

It satisfies |^ = 1, and the third component is 

PiK) = Vi-kP. (16) 

This is real valued for |it| < 1, meaning that the waves are propagating and not evanescent. In general there are forward 
and backward propagating waves, as explained in Section 4.4, but in our scattering regime the waves move forward, 
which is why we have P{k) with positive sign in (15). 

The vectors 

u{k) = (yS(*:)j^, -\k\), a^{K) = (j^,0), (17) 

where - {-K 2 ,ki), distinguish the two modes associated with wave vector k. Note that the triplet {u{k),u^{k),k\ 
is an orthonormal basis of Thus, equations (13) and (14) are decompositions in transverse waves, which are 
orthogonal to the direction of propagation along £ The amplitudes a(K,z) are for the transverse magnetic modes, 
because they do not contribute to the longitudinal component of the magnetic field. They multiply the vector u'‘~(k) 
in equation (95). Similarly, a-‘-(K,z) are the amplitudes of transverse electric modes, because they do not contribute to 
the longitudinal electric field. 

For any wave vector il, the electric and magnetic plane waves 

£(/f,z) = p^^(K)[a(K,z)it(K) + a-^(K,z)ii'^{K)], 
lK(/f,z) = ^^'^P^HK)[a{K,z)ii'^iic) - a-^iK,z)iiiK)], 

are orthogonal to ic and to each other. The amplimdes a{K, z) and a-‘-(/c, z) are random, and encode the scattering effects 
in the medium. In the system of coordinates with axes along u(/c), u^{k), and k, which is right-handed because 

U{k) X U^{k) - li, iix UiK) - X i? = !?(*:), 


we can quantify the evolution of the state of polarization of the waves using the coherence matrix 

W6i(/!:,z) -H 62 (#t,z) ©3(*:,z) -h ! 64 (*:,z)\ 

2 \ 63 (*:,z) - i&A(,K-,z) 6 i(#i:,z) - &2{k-,z))" 

where t denotes complex conjugate and transpose, and 

6i(*:,z) = E[|a(#t,z)|^ -r \a^(K,z)?-], &2 {k,z} = E[|fl(A:,z)p - \a^iK,z)?'], 

& 2 (i^,z) - 2Re(E[a(#t,z)a-^(/!:,z)]), 64 (*:,z) = 2Im(E[fl(*:,z)a-^(A-,z)]), 


J’(/f, z) = E 


a{K, z) 
i^{k,z) 


a{K, z) 
2^{k, z) 


are the components of the Stokes vector &{k,z) - (©i(/!:,z), &2{k,z), 63 (*:,z), © 4 (iit,z)). We use henceforth the bar to 
denote complex conjugate. The degree of polarization is defined by 


Pol(#i:, z) = 


{& 2 (K,zf + & 2 {K,zf + ©4(^,Z)")'^" 


<Sl{K,z) 

It is a number in the interval [0,1], due to the Cauchy-Schwarz inequality 

|2 


|E[a(*:,z)a-^(A',z)]| < E[|fl(K:,z)P]E[|fl'^(A-,z)P]. 


(19) 


The wave is fully polarized when Pol(#t, z) = 1 and unpolarized when Pol(*:, z) = 0. The mode amplitudes of fully 
polarized waves have a deterministic linear dependence so that 

E[|a(Af,z)|^]E[|a-^(/c,z)|^] = |E[a(A-,z)a-^(#i:,z)]| . 

The mode amplitudes of unpolarized waves are uncorrelated and have the same mean power. 
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2.4. Main results 

If the medium were homogeneous, the modes would be independent, with constant amplitudes Goiic) and a^iK) 
dehned by the current source (10), 


ao{K) = 




/X\ 

dxJz[-)e 


atix) = 




f 

Jk2 

n- f 

kl Jr 2 


K 


-■ f 

Jr2 


dx 


-ikK-z 


-ikK-i 


( 20 ) 

( 21 ) 


In the random medium, in the scattering regime dehned by (11-12), the mode amplitudes {a{K, z), a-‘-(/(, z)) for |iit| < 1 
and z L form a Markov process whose statistical moments can be characterized explicitly. We are interested in the 
hrst two moments which describe the loss of coherence of the waves and the transport of energy by the modes. 


2.4.1. The mean mode amplitudes 

Let us denote the vector of the expectations of the mode amplitudes at wave vector k by 


A{k,z) - 


I E[a(*:,z)] \ 
\E[a-^(*:,z)]/' 


Its initial value in the source plane z = 0 is 

.4(..o) = AW = (»;«), 


( 22 ) 


(23) 


and for z > 0 we have 

A(k, z) = exp [Q(/i-)z] Ao(/!:). (24) 

Thus, the scattering effects in the random medium do not average out, and the range evolution of the mean mode 
amplitudes is determined by the matrix Q(#t) in the exponential. This complex symmetric matrix is given by 


Q(^) = • 




f 


d(kK') 

{Inf 


r(*:, k)T{k‘ 


',k) f dx f 

Jk2 Jo 


c/zlk 


2 0/’ 


for / dehned as in (15), £ - (x, z), and 


r(^,^') 


' kik'l 

a/swaT) 



w'\ V m 


_£ ^ I Pin) ' 

k| ■ k'l -y /?(*:') 

_ 1 _ 

1*1 l*'l sJI3(k')I3(k), 


(25) 


(26) 


Note that r(iit', k) = r(/i’, K'f and since the autocorrelation 3? is an even function, and the power spectral density IR is 
non-negative, the real part of Q(*:) is a symmetric, negative dehnite matrix 


Re[Q(^)] 




f 

Jk'Kl 


d(kK') 

~(2iff 


r(A:, K')r'^{K, k') Ji{ki{i^ - ic)). 


(27) 


In the case of transverse isotropic statistics of the media, where lR(f^ with P = (r, rp depends only on |r| and r^, 
we obtain by integrating in (25) over k' in polar coordinates that Q(*:) is diagonal and moreover, it only depends on 
litl. Therefore, the mean wave amplitudes decouple and they decay exponentially in z on the range scales (scattering 
mean free paths) 


*S(^) = -l/Re[Q„(^)], ^^(^) = -l/Re[Q 22 (^)], (28) 

which depend only on |*:|. We refer to Section 8 for illustrations with numerical simulations in statistically isotropic 
media, and to Section 9 for the analysis of the high-frequency limit and connection to the white-noise paraxial regime 
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described in Appendix C. We also show in Appendix B that, if is isotropic i.e., only depends on then 
Re(Q(/i:)) is proportional to the identity matrix, so the scattering mean free paths (28) are equal. 

In general media the mean mode amplitudes are coupled, but they still decay exponentially in z. Indeed, matrix 

S(*:) = -Q(*t) - Q(a:)^ = kjEa f T{k,k')T'^ it)) (29) 

4 (27r)2 

is real symmetric and positive definite. Therefore, the mean amplitudes satisfy 

d^\\A{K, z)\f' - A{k, z)"*' [Q(#t)''' + Q(/!:)jA( k:, z) = -A{k, z)'''S(#t)A(*:, z), (30) 

and if we let Ai(a:) > A. 2 {k) > 0 be the eigenvalues of S(a:) and use Gronwall’s lemma, we get 

The scales 2 /Ai(a') and 2/A2(it), on which the mean wave amplitudes decay exponentially in range, are the scattering 
mean free paths. 

The exponential decay of the mean mode amplitudes, and therefore of the mean (coherent) fields E[.E(®)] and 
E[i3'(®)], is a manifestation of the randomization of the waves due to scattering. Energy is conserved, as shown in 
Section 5.3, 

r z)l^ + |fl'^(/i:,z)l^l = constant, (32) 

Jw<i ^ J 

so the second moments cannot all decay. Therefore, the mode amplitudes have large random fluctuations. 

Remark 1. The exponential decay of the coherent field is a general phenomenon in random media, and can be 
characterized under many scattering regimes, including when the forward scattering approximation is not valid. We 
refer to [18] for an overview. 


2.4.2. The transport equations 

The energy density (Wigner transform) of the modes at wave vector k is defined by 


'W(*:, £) — 


f 


exp [ikq ■ (x + V/3(q)z)] E 


a{K+ |,z) \l a{K- |,z) 


a^{K+ ^,z) \a-^{ic- |,z) 


a? = (a;,z). 


(33) 


where the integral is over all q such that |a: + q/2\ < 1. It satisfies the transport equation 

d^V^iK, £) - VjS(/!:) ■ Va;'W(iit, £) — Q(iit)'W(A', £) + 'W(*:, a()Q(*:)^ 

f y)W(y , £)r(/t', #t) iuiK -k)], 

Jk'Ki (221)2 V J 




(34) 


for z > 0. The derivation of this equation, given in Section 7, is one of the main results of this paper. We show 
some numerical illustrations in Section 8 in statistically isotropic media and relate it to the radiative transport theory 
in [13, 9, 10] in Appendix B. We consider the high frequency limit Ajt ^ 0 in Section 9, where we can connect the 
results to those in the white-noise paraxial regime [8] as described in Appendix C. 

The integral over x of the Wigner transform is the Hermitian, positive definite coherence matrix 


‘fiK, z) = E 


a{K, z) \ / a(K, z) 
a^{K,z)]\a^{K,z) 


f 

Jr 2 


dx 'W(*:, £). 


(35) 


Its diagonal entries are the mean mode powers, and the relation with the Stokes vector is as in (18). The coherence 
matrix evolves from its initial value TiK, 0) = A(k, 0)A{k, 0) ’ at the source according to equation 

kp' p ftt kid 1 — / \ 

dyP{K,z) ^Q{k)9{k,z) + ‘P{k,z)Q{k)' + —-— I ^ T(k,k')9(k',z)T{k',K)‘T i[k{(il- it)] . (36) 

4 {2ny '' 2 
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It follows by direct calculation, using the commutation properties of the trace and (29), that 




( 37 ) 


which is consistent with the conservation identity (32). Moreover, (36) can be written in integral form as 




k^e^c 


f */ 

Jo J|a :'|<1 


d(kjc') 

'(W 


•R[k{(K-K))x 


x[eQWfe-^')r(^^ it')]?*)*:', z')[eQ('')("""')r(A:, a:)] ', 


where we used the symmetry relation r(A:',A') = r(A:, it')^- Since J’(a:, 0) is Hermitian positive definite by definition, 
and the power spectral density Jl is non-negative, we conclude that equation (36) has a Hermitian positive definite 
solution for all z, which is consistent with (35). 


3. Mathematical formulation of the problem 

In this section we give the mathematical formulation for the analysis of the electromagnetic wave field in the 
random medium with electric permittivity (5) . We begin in Section 3.1 with the derivation of the 4x4 system of 
partial differential equations satisfied by the transverse components of the electric and magnetic field. Our goal is to 
analyze the solution of this system in the scattering regime (11-12), made precise in Section 3.2. This scaling defines 
the asymptotic regime in which we can characterize the statistics of the mode amplitudes using a Markov limit. 


3.7. The equations for the transverse waves 

Let us eliminate from Maxwell’s equations (1) and (2) the longitudinal components and H^of E - {E,Ef) 
and H - (H, Hf), and derive a closed system of partial differential equations for the transverse electric and magnetic 
fields E and H. We obtain that 


HfdS) =-V-" ■ E(£), 

(38) 

tO/Jo 


EA^) = ■ Hitf) - 4 • jm. 

(39) 


where is the unit vector along the range axis, V = (dx,,dx 2 ) is the gradient in the transverse coordinates and 
V-*- = {-dx 2 ,dxt) is its rotation by 90 degrees, counterclockwise. The transverse components satisfy 




JfxlX) 


_ 

'' X 

(Jj 

e(£) 


d^E{x) - ■ 


= -iioei^)E(cg) -V-"[V-" ■ E(£)] + J{xlX)5{z), 

(Dfio 




(40) 

(41) 


with - (-7/2,77i), the rotation of iT = {H\,H 2 ) by 90 degrees, counterclockwise. We study the 4x4 system 
(40-41), written in more convenient form in terms of the scaled, rotated magnetic field 


U{£} = 


(42) 


The scaling by the impedance ensures that E and U have the same physical units. Equations (40-41) become 


d2,Ei£) = ikU{£) -H -V 
k 



Jzix/X) 

6{z), 

(43) 

Sr(£) 



d,U(cg) = iksf£)E(£) -H [V-" 

£;(al)] 

- J(xlX)6{z). 

(44) 
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3.2. Asymptotic regime 

We model the separation of scales in (11-12) with two dimensionless parameters e and y, ordered as 

0 < e « y < 1, (45) 

and defined by 

A A 

“V 

Our primary asymptotic analysis is for e ^ 0. The parameter y remains finite, except in Section 9 where we study 
the high-frequency limit y —» 0 in which the transport equations simplify and become equivalent to those in the 
white-noise paraxial regime. 

We let L be the reference length scale, and introduce the scaled length variables 

x'^xlieL), z^zlL, L'^LIU X'^XlL^ely,. (47) 

The scaled wavenumber is k' - kLe — 2n and the normalized standard deviation of the fluctuations is 

a' = = 0(1). (48) 


As we will see below, the dimensionless parameter 


jj = ytjX = 0(y), 


(49) 


controls the opening angle of the cone (beam) emitted by the source. 

Substituting (47) into (43-44) and dropping all the primes, since all the variables are scaled henceforth, we obtain 


dzE%x,z) 


ik 

1 . j 

ia 

r / 7\ , 1 

e 

I+^V(V.) 


v^ya;,y-j V ■ U {x,z) 



' [yx,y-^X ■ U'^{x,z) 


i 

k 


'^xJzi7jX)6iz), 


(50) 


and 


dzU^(x,z) 


ik 

€ 


i + ^v-(v-.) 


ikcY I .z \ 

E^(x, z) + —qy^v\yx, y-j E^(x, z) - J{yjX)6(z), 


(51) 


for 0 < z < L and I the 2x2 identity matrix. For ranges z < 0 and z > L the equations are simpler, as all the terms 
involving the process v vanish. The fields and C/'^ are approximations of E and [/ for e « 1, and equations 
(50-51) follow by neglecting terms of order and higher in (43-44). 


4. Wave decomposition 

Because the interaction of the waves with the random medium depends on the direction of propagation, we de¬ 
compose E^ and over plane waves, using the Fourier transform with respect to the transverse coordinates x. We 
denote this transform with a hat, to distinguish it from the Fourier transform with respect to the three dimensional 
vector A - (x,z). We have 

E\kK,z)^ f dxE^(x,z)e-‘^‘''^, (52) 

Jr 2 

and similar for U^{kK, z), where k is the scaled wave vector. The inverse transform is 

E^(x,z)^ f ^E%kK,z)e‘'‘'^-^. (53) 

Jr 2 (^ttY 

We begin in Section 4.1 with the decomposition in homogeneous media and then consider random media in Section 
4.2. We state the energy conservation in Section 4.3 and interpret the wave decomposition in terms of transverse 
electric and magnetic plane waves in Section 4.4. 
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4.1. Homogeneous media 

The transformed fields in homogeneous media are Eg{kK,z) and Ug(kK,z). They satisfy the system of ordinary 
differential equations 


E^(kK,z) 

[U^(kK,z)j 


^ ik 
= -M(it) 


E^ikK,z) 

[U‘(kK,z)j 


S(z) 




0 


I-/tc 


with 4x4 matrix 

" (l - 0 

and symbol ® denoting vector outer product. The Fourier transform of the current source S = (J, J^) is 


(54) 


(55) 


J(kK) 


-f 

Jr 2 


dx J{x)e 


~ikK-x 


JzikK) 


-f 

Jh 2 


dx Jz{x)e 


~ikK-x 


(56) 


Naturally, the fields are outgoing and bounded away from the source. Moreover, the Fourier transform ( J(kK), J^ikK)) 
of the current source is supported at |iit| < 1 /k. Therefore, the excitation in equations (54) is supported at transverse 
wave vectors satisfying 

kl < T < 1 - (57) 

k 

The plane wave decomposition is based on the diagonalization of the matrix M(iit). This has two double eigenval¬ 
ues ±/1{k) and the eigenvectors are 







and 




(58) 


They are linearly independent when |a:| 4 1, so they form a basis of in which we can expand the solution of (54), 


'K(kK,zi 

JU^{kK,z), 


Vo{K,z)f+(K) + V^{K,z)f+(K) +Xo(K;Z)lp-{K) + xt{K,z)lf/t{K). 


(59) 


Because the frequency is fixed, we do not include the wavenumber k in the argument of the coefficients in the expan¬ 
sion. Equation (59) is the Flelmholtz decomposition of the electric field and the rotated magnetic field, because the 
components along k correspond to curl free vector fields in the x domain, and the components along to divergence 
free fields. 

The expressions of the coefficients in (59) are obtained by substituting into (54) and using the linear independence 
of the eigenvectors. We obtain that 

Vo{k,z) = Xo{K;Z) = fio(#t)e^ 

and similar for Vg and;^^^. Consequently, the electric field is given by 

E^g(kK,z) = l(o,co)(z)[fl„(^)y6'^'(^)^ 

1^1 1^1 ■* 

+ l(-oo,0)(z)[ - (60) 

1^1 1^1 ■* 

and the rotated magnetic field is 

Ug(kK,z) = l(0,oc)(z)[ao(A:)/3^'^^(/f)-^ -t- 

|/t| |/t| -* 

|A!| |/t| -* 
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( 61 ) 












Recalling the Fourier transform (52), we see that this is a plane wave decomposition with wave vectors {k, ±/3(k)) 
multiplying {x,zle) in the phases of the exponentials. The plus sign corresponds to forward going waves, along the 
positive range axis, and the negative sign to backward going waves. 

In (60) we used the radiation conditions, so that the waves are outgoing and bounded, and the amplitudes are 
determined by the jump conditions at the source 

._. ._. _ _. / \ ._. 1 ._. / \ 

ElikK, 0+) - ElikK, 0-) = - , U^„{kK, 0+) - U^,{kK, 0-) = -^ J - . 

jj \yjl Jj \7 jI 


We obtain that 


aoiic) = 


bo(K) = Z-T 

2 rf 





Ir.,/ 




{Jj! 


?l/2/ 


?l/ 2 / 


kK 


kK 




(62) 

(63) 


Note that equations (62) written in dimensional units are the same as (20)-(21). Moreover, as noted above, the 
excitation in equations (54) is supported at transverse wave vectors satisfying (57). Thus, there are no evanescent 
waves in the decompositions (60-61). 


4.2. Random media 

The Fourier transforms of the wave fields in the random medium satisfy the system of equations 


dz 

'E^(kK,z) 

ik 

= -M(/t) 

{E^(kK,z)\ , 

Cz +1(0 L)(z) 


'E'^\ 


JJ^(kK,z), 

e 





{kK,z), 


(64) 


derived from (50-51), with radiation conditions at z < 0 and z> L, and source conditions at z = 0. The leading order 
term in the right hand side involves the same matrix M(iit) as in the homogeneous medium. The perturbation due to 
the random medium is given by the integral operator defined by 


AT 





(k/t, z) 


ika r d(kK') k(/c -/c') 

gi/2y2 J (27r)2 \ y 

ika^ r dikid) -:,(k{K- k') z\/0 

y2 J (27r)2 \ y ’^e)\0 


E^{kK',zS 
0 ][uHkK',z), 


E^ikid,zS 
0 ][uHkK',z), 


(65) 


where 7 and are the Fourier transforms of v and with respect to the first argument. 

We use a similar decomposition of the waves as in the homogeneous medium, based on the same eigenvector basis 
that diagonalizes the matrix M(it) in the leading term of (64). We obtain that 


'E‘(kK,zi 

JJ^ikK,z), 


-ya\K,z)il/+(ic) + a^’-^(*:,z)iAi:('t)]e'^***^ + \b\K,z)^SK) + (k, z)>j/i(K)^e 


( 66 ) 


where the amplitudes a^iK,z), a’^'^{K,z), E{k,z), and b^’-'-{K,z) are no longer constants determined by the current 
source density, but random fields. We keep all the components of the waves, forward and backward going, because of 
scattering in the random medium. 

The radiation conditions are 


a^{K,z) - a‘^’-^{K,z) - 0, if z < 0, and E(k,z) - b^’'^{j(,z) - Q, if z > T, (67) 

because the medium is homogeneous outside the range interval (0, L). This also implies that 

a'^{K,z) = a^{K,L), and a^’^{K,z) - a^'^{K,L), ifz>T, ( 68 ) 
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and 


b'^{K,z) - b^{ii,0-) and b'^’^(K,z) = b‘^'^{K,0-), ifz<0. 
The excitation comes from the jump conditions at the source 


(69) 


E^{kK, 0+) - E^{kK, 0-) = -r/J — 

jj \yj 


. _. ._. 1 ._./ 

U%kK, 0+) - U%kK, 0-) = —r J — 

7j Uj 


which give the following initial conditions at z = 0, 

0 +) = ao(ic), a'^’-^(K,Q+) = a^iK), 

and 


(70) 


E(k, 0-) = ba(K) + E(k, 0+), 0-) = b^(K) + E'^{k, 0+). (71) 

Equations (70) say that the forward going waves leaving the source are the same as in the homogeneous medium. 
This is physical, because the waves must travel a long distance before they are affected by the small fluctuations of 
the electric permittivity. Equations (71) say that the waves at z < 0 are given by the superposition of those emitted by 
the source, modeled by boiK) and bg(K), and the waves backscattered by the random medium, modeled by E(k, 0+) 
and Z7^ -^(/f,0+). 

To determine the amplitudes in the random medium, we substitute equations (66) into (64), and use the linear 
independence of the eigenvectors and If we let 


Y%k,z) 


' a^(K, z) ' 
a‘^’-^(K,z) 
b%K,z) 


(72) 


we obtain that 


dY^ 

dz 


{k,z) = 


ika 

2'y2g-l/2 

ika^ 


/ 


d(kK')^lkiK- k') z 


I 


(271)2 
d(k/(') — 


2y2 J (27r)2 


r 

k(A: - k') 
7 


>r- 


-.7- 


¥(k,i^,-^Y%k',z) 

;(«./, i)FV, A 


(73) 


in z e (0,L), with boundary conditions (67) and (70). We are interested in the propagating waves, corresponding to 
|*:| < 1 in (73), and we explain in Section 5 that in our regime the evanescent waves may be neglected. This means 
that we can restrict the integral in (73) to vectors id satisfying l/t'l < 1. 

The mode amplitudes are coupled in equations (73) by the 4x4 complex matrices Fj/t, k' ,() and G)*:, k' with 
block structure 


nK,K\o 




(74) 


and 


G)*:, k ',0 




(75) 


where the superscripts of the 2x2 blocks indicate which types of waves they couple. The dependence on ^ of the 
blocks in F can be factored as 


F'“=(*:, , 0 = r'“=(*:, ( 75 ) 

F'’\k, k',0 = ^ ( 77 ) 

F‘‘\k, k',0 = Y\ic, *;')e-'4-8('^')+/sw]f ^ (78) 

F*"(7i:, k',0 = (79) 
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with 2x2 real-valued matrices 


and 


T^‘'{k,k) = 


T‘’\k,k') = 


kik'l 






1*1 I*'I -y /?(*) 






w |*'| >//?(*') 

_ 1 _ 

w I*'I y/l3(K')0(K). 


_* ^ / jSW 

1*1 ■ i*'i -y /?(*') 

K k' 1 


1*1 I*'| -y /J(*) 1*1 I*'| 

on the diagonal. Note that is the same as the matrix F in (26). The off-diagonal matrices are 


and 


F^^a:,*:') = 


F*"(x, x') = 


yjl3(K)l3(K') 1*1 I*'I 


1*1 I*'I -y /?(*) 


1*1 |*'IA//?(*') 

_ 1 _ 

W I*'I yll3(K')l3(K). 


kik'l 


- + 1*1 ■ I*,I ylf^(KWK') 1*1 ■ *^,| 


-y,8(*),8(*') 1*1 I*'I 


1*1 I*'I -y /?(*) 


The diagonal blocks in G are 


G‘“‘{k,k',0 = 


kik'l 


V/3(*)/3(*') ^ 

0 0 


1*1 k'ly/?(*') 

_1_ 

1*1 I*'I V;3(*')j8(*) 






kik'l 


and the off-diagonal blocks are similar 

G‘‘^iK,K',0 = 

G''"(a:,a:',^) = -G“*(x,x',0- 


V<8(*)A*') 

0 Oy 


-ii:[/S(*')+/t(*)]f 


(80) 


( 81 ) 


(82) 

(83) 

(84) 

(85) 

( 86 ) 

(87) 

( 88 ) 

(89) 

(90) 


4.3. Energy conservation 

The diagonal blocks of coupling matrices F and G satisfy the symmetry relations 
F‘’“(a:, k', 0^ F^^iK', K, 0\ G‘“‘(k, k’,0 = G“‘'(a:', k, 

F‘’\k, k',0^ F^’^k, k, G‘’\k, k,0 = G^'^i/, K, 

Using these in (73) we obtain the conservation identity 

r - \b''(ic,z)f- - \b'^’-^(K,z)f] = constant, 

J|*l<i (27r)^ L j 

independent of z. 

Now let us recall that ^Ke[E'^ x iT^] is the time average of the Poynting vector of the time harmonic electromag¬ 
netic wave [19]. Its z-component is 

--Re[£'^(a;,z) ■ H’'-^(x,z)\ - -|-Re[£'^(a:,z) ■ 17^(a;,z)j, 

where we used (42). The integral of this is the energy flux which must be conserved 

(7a:Re|^£l^(£c,z) ■ I7'^(a:,z)j = constant. (91) 


Jr2 


Our identity (90) is the statement of this energy conservation, written in terms of the wave amplitudes. It follows from 
(91) by substituting equations (66) and using PlanchereTs theorem. 
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4.4. Transverse electric and transverse magnetic waves 

To interpret the wave decomposition (66) as an expansion of the electromagnetic wave field in transverse electric 
and magnetic waves, we recall the definitions (38)-(39) of the longitudinal electric and magnetic field. After the 
scaling, and in the Fourier domain, these become 


H^,{kK,z) = C ■ E%kK,z), 

El(kK,z) - ■ H^{kK,z) - -K ■ U%kK,z). 


(92) 

(93) 


Then, using (66), the definition (58) of the eigenfunctions, and the inverse Fourier transform (53), we obtain that the 
electric field E^) is given by 






#(-,z)= r HK)i[a'^(K,z)it(K) + a''’^(K,z)it'^(K)]e' 

e Jw<i (27r)2 ' 

+ [ - b%K, z) u~(k) + z) u^{K)\e' 

and the magnetic field is 

e Jm<i (27ry ' 

+ [E(/c,z)u'^(/() + E"^(/c,z)u~(/()]e‘^^'‘ 


(94) 




(95) 


for £ = (x,z). Here we used the wave vectors ^ = {k,I3{k)) and ic - {k, -j3{K)) of the forward and backward going 
waves, and 

«(*:) = «“(*:) = h-^(/!:) = ^ oj . 

The triplets {u{k),u^{k),k\ and {u~{k),u^{k),k ) are orthonormal bases of K^, when |*!:| < 1. Thus, (94) and (95) are 
decompositions in transverse waves, which are orthogonal to the direction of propagation along the wave vectors ic 
and iT, respectively. The amplitudes a^iK,z) and b^{K,z) are for the forward and backward transverse magnetic plane 
waves, and a^’-'-iK,z) and b'^’-‘-(K,z) are the amplitudes of transverse electric plane waves. Note that the forward going 
part, written in dimensional units, is the same as in (13) and (14). 


5. The Markov limit 

Here we describe the e —> 0 limit of the random mode amplitudes satisfying the two-point linear boundary 
value problem (73) with boundary conditions (67)-(70). We begin in Section 5.1 by writing the solution in terms of 
propagator matrices. Then we explain in Section 5.2 why we can neglect the backward and evanescent waves. The 
limit of the forward going wave amplitudes is given in Section 5.3. 


5.7. Propagator matrices 

The 4x4 propagator matrices z; Kg) are solutions of the initial value problem 


dP^{K,Z-,K„) 

dz 


ika 

2y^eF2 

ika^ 

2y^ 


f 

Jk'Kl 

f 


d(kK') k(K - k') z 


-.7 

r e 


d(kK') / k{K - id) z 

7^’^ \ r 


f(^,^',-)pV,z;^«) 

- ,z-,Ko), 


(96) 


with initial condition P'^(*:,z = 0; Kg) - 6 {k - Ko)\, where I is the 4x4 identity matrix. The solution of (73) satisfies 

Y%k,z)^ f diCgP^(ic,z;iCg)Y^(iCg,0), 

JkJ<i 
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for any z e [0, L\. In particular, at z = L and using boundary conditions (67)-(70), 


' a'^(K,L) ' 
L) 

0 

0 j 


ao(,Ko) \ 



dKo P‘^(a-, L; Kg) 


atiKo) 

b\K„,0) 


(97) 


5.2. The forward scattering approximation 

The e —» 0 limit P of P'^ can be obtained and identified as a Markov process, meaning that it satisfies a system of 
stochastic differential equations. We refer to [11, 12] and Appendix A for details, and state here the result. 

In the limit e ^ 0, the terms of order one in (96) i.e., those involving the kernel G, become equal to their average 
with respect to the distribution of v and to the rapid 0(1/e) phase. We denote the average with respect to this phase 
by {■), and obtain 


ika^ f d(kK') _ f':; / k(K -id) 

2r Jk'i<i (2^)^ [ \ r 


iJcct 

{G(k,k',-))P(k',z;Ko) = ^-D?(0)<G(/!:,/f,-)>P('t,z;>fo)- 


Due to the large phases of the off-diagonal blocks in G we have 


<G(#t, K, •)) 


(-1 

kP 0 

m 0 
lo 


0 0 Of 
0 0 0 
0 1 0 
0 0 0 , 


(98) 


a diagonal matrix, so the order one terms in (96) do not introduce any wave coupling. 

The order C>(e '^^) part in (96) i.e., involving F, gives diffusive terms. Let us split the propagator matrix into four 
blocks: 

p6ffl,6 pfl/ 7,6 
pi7fl,6 'pbh,e 

and consider first only the propagating waves. The stochastic differential equations for the limit entries of z; Kg) 

and P*‘’’^(a:, z; k„) are coupled to the limit entries of P‘^‘‘’‘^(k, z; Kg) and P**’'^(/!:, z; Kg) through the coefficients 

Jk \ r / Jb 3 r 




Here fR is the Fourier transform of fk with respect to the cross-range variable, and fk is the three-dimensional Fourier 
transform (8). The sum j3(k) + /?(*:') in the phase of (99) is because the phase factors in the matrices F"* and F*" 
are ±(fi(K) - 1 - f5(K'))f. The coupling between the entries of P“‘’’^(*:', z; k„) and between the entries of P^’*’'^(it', z; Kg) is 
through the coefhcients 

r t/r; ik ( ~ ^ \ jr] ^ f = 4 

Jb \ r / Jk3 


k(K - k') 
y 


( 100 ) 


because the phase factors in matrices and F** are ±(J3(k) - f5(K'))f. 

We conclude that the interaction of the forward and backward going waves depends on the decay of the power 
spectral density. Recall that the source gives wave amplitudes supported at k| < < 1 and fk is negligible outside 

the ball with center at 0 and radius one in K^. From (49) and (57) it is possible to choose some e (jj/k, 1) such 
that y satisfies 


r 


(101) 


Under these circumstances, for all k' satisfying k'l < k^, the coupling coefficients (99) vanish, because 



r 


- [k - + Ifi(K) +f(K')ff^^ > - \P(K) +P(k')] > 


2kA^M) 

r 


> 1. 
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This means that the forward and backward wave amplitudes are asymptotically decoupled as long as the energy of the 
wave is supported at the transverse wave vectors k satisfying l/tl < k^. 

Nevertheless, the forward going amplitudes are coupled with each other, because 

k\K-K'\ kv o '.nl/2 

J -1 = - k - + {J3(K)-li{K')f\ 

■y '■y L J 

is small for at least a subset of transverse wave vectors satisfying |iit|, l/t'l < a:„. Due to this coupling there is diffusion 
of energy from the waves emitted by the source with l/tl < jjlk, to waves at larger values of |iit|, as explained in 
Section 8. This is why we take km > jjIk in (101). By assuming that a^(K,z) and are supported at |*:| < km 

we essentially restrict z by Z^, so that the energy does not diffuse to waves with |*:| > a:„ for z < Z^^. Physically, it 
means that the three-dimensional wave vectors k of the forward going waves remain within a cone with opening angle 
smaller than 180 degrees. We will see that the evolution of the *:-distribution of the wave energy is described by a 
radiative transfer equation, which means that the wave energy undergoes a random walk (or diffusion). We will also 
see that the diffusion coefficient is of the order ya^. Thus, the *:-distribution of the wave energy reaches after a 
propagation distance of the order of such that a^yZ^ - . Recalling the scaling (48) of the standard deviation of 

the fluctuations and (12), we see that it is possible to choose Z^ = Oil) (i.e., similar to L in dimensional units). 

The evanescent waves can only couple with the propagating waves with wave vectors of magnitude close to 1. 
Thus, as long as the energy of the wave is supported at wave vectors satisfying l/tl < km, assumption (101) implies 
that the evanescent waves do not get excited. 

The forward scattering approximation amounts to neglecting all the backward and evanescent waves in the system 
of stochastic differential equations. It is justified in the limit e ^ 0 by the decoupling of the waves, the zero initial 
condition of the evanescent waves, and the zero condition at z = L of the left going wave amplitudes. 


5.3. Markov limit of the forward going mode amplitudes 

Under the forward scattering approximation, the wave amplitudes a'^(K,z) and a'^’-‘-(K,z) satisfy the initial value 
problem 


/ a^(K,z) \ 

dzW’^(>c,z)j 


ika 

2'y2^1/2 

ika^ 

2y^ 


f 

f 

Jk'Kl 


diicK'), 
d(kK') ■ 


k{K - k') z 

->r- 

r e 

k{K -id) 


7 


,7-\G^ 




a^K',z) 

2^’^(k',z) 

a^(K',z) 

2^'^{k',z) 


( 102 ) 


for z > 0, and the initial condition 


/fl"(^,0)\ 

K'^(^,0)j 


Ao{k), 


(103) 


with Ao defined in (23). These equations conserve energy, as the special structure of the kernels F'*" and described 
in equations (76), (84), and (88) implies that for all e > 0, for all z > 0, 



d(kK) 


[|a^(*:,z)|^ 


Jw<l 


d(kK) 


[|ao('t)l^ 


1^0 ('t)l^] ■ 


(104) 


We refer to Appendix A for the details on the Markov limit e —> 0. In particular, we explain there that the process 




' Re(a'^(#i:, z))' 
Im(a'^(#t, z)) 
Re(fl^’-"(it,z) 
Tm(a'^’-^(#i:, z). 


KeO 


for 0 = {/t e K^, |a-| < 1), 


converges weakly in C([0, L],D') to a Markov process X(z), where £)' is the space of distributions, dual to the space 
£)(0,K^) of infinitely differentiable vector valued functions in K^, with compact support. The generator of X(z) is 
given in Appendix A, and we denote henceforth the limit amplitudes by a(K,z) and a-‘-(*:,z). Their first and second 
moments are described in the next two sections and are summarized in Section 2.4. 
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6. The coherent held 


The vector 


A{k, z) - 


/ E[a(*:,z)] \ 


= limE 


/ a'^(/(,z) \ 


(105) 


whose components are the expectations of the limit amplitudes of the modes at wave vector satisfies as explained 
in Appendix A the initial value problem 


d:A(K, z) - Q(ic)A(k, z), z> 0, 

with initial condition A(k, 0) = Ao(k), and 2x2 complex matrix Q(*:) given by 


Q(^) - ■ 


4^3 


r d{kK') , . 

—^T{K,K)TiK,K) d^'R 


k{K -id) 

r 


,ae 


(1 O' 

2 0, 


In dimensional units this is the same as (25), with F given by (26), which equals T®" defined in (80). 
As stated in Section 2.4.1, the real part of Q(a:) is a symmetric, negative definite matrix 


2^,2 


Re[Q(^)] = 


k^a- 

8y^ 


r dikK’) ~ 

'k(K-K'y 

i|.'|<i (2;r)2 

[ 7 ) 


r(ic, K')r(K',K), 


(106) 


(107) 


(108) 


because Hi is non-negative and r(/f',#t) = r(*:, a:')^. It is an effective diffusion term in (106), which removes energy 
from the mean field and gives it to the incoherent fluctuations. This causes the randomization or loss of coherence of 
the waves. The imaginary part of Q(a:) is the sum of two terms. The first, given by 


kV r 

4r^ Jk'Ki 


d{kK') 


r 




7 


r(A:, A■')^(A■^ *■), 


is an effective dispersion term, which does not remove energy from the mean field and ensures causality'. The second 
term in (107) is a homogenization effect due to the part in (102). It accounts for the slightly different velocities of 
the transverse electric and transverse magnetic waves, with discrepancy of the order of 0(e). 


7. The transport equation 

To describe the transport of energy, we study the z-evolution of the Hermitian positive definite coherence matrix 


(P(k,z) = limE 


a^(K, z) W a'^iK, z) 
a'^'-^(K,z)j\a'^’-^(K,z) 


(109) 


The diagonal entries of (P are the mode powers. The evolution equation is derived as explained in Appendix A, 


d^(P(ic, z) -Q(ic)(P(k, z) + RiK, z)Q(a:)''' h- 


k^a^ 

4y3 


/ 


d(kK') 

~(2^ 


r(A-, k')(P(k' , z)T(k' , k)R. 


k(^-K')'^ 
7 


( 110 ) 


for z > 0, with initial condition ?*(/<•, 0) = Ao(ic)Ao(icp■ In dimensional units this is the same as (36). 
We have from (104) that 

f ^ f p^\\ao(id)?'+ \at{id)?'\ 

Jk|<i (2^)^ J Jk|<i (27r)^ L J 


^ If we write the coherent wave fields in the time domain, using the inverse Fourier transform with respect to the frequency to, we obtain a 
causal result. 
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for all z > 0 and 0 < e « 1, so energy conservation must hold in the limit e ^ 0. This is derived from (110) using the 
commutation properties of the trace, as stated in (37). 

By studying the second moments of the wave amplitudes, we can also quantify the decorrelation of the waves due 
to scattering in the random medium. The calculation of these moments uses the generator in Appendix A and the 
result is that the wave amplitudes at different transverse wave vectors k + k' sit decorrelated 


limE 


/ a^K,z) W a'^(K',z) \ 

\a<^-^(K,z)j \a'^’-^(K',z)j 


A(k, z)A(k', z)^ ■ 


Intuitively, this is because these waves travel in different directions and see a different region of the random medium. It 
is only when \k- k'\ - 0(e) that the waves are correlated, and we can calculate the energy density (Wigner transform) 


W(*:, x,z) 


=iim r 

J 


d(kq) 


exp 


^ikq ■ (Vy6(#t)z + a:)jE 


a^{K+ ^,z)\( a^(K- ^,z) 


a'^’-^{K+ ^,z))W’-^(k- ^,z) 


( 111 ) 


It satisfies the transport equation 

d^WiK, x,z) - ^/i(K) ■ V3 ,W(a:, x,z) - Q(k)W(k, x,z) + 'W(*:, a;,z)Q(A:) 

k^a^ r d(kK') , , , ~(k(il-K‘ 

4r J\K'\<1 (27r)2 y 


( 112 ) 


forz > 0. In dimensional units this is the same as (34). When the initial condition is smooth, we have from (111) that 


W(*:, X, 0) = 6(x)Ao{k)Ao(k)\ (113) 

and therefore at z > 0 

W(/t, x,z) - 6(x + V/3(a:)z) (P(ic, z), (114) 

so energy is transported on the characteristic 

K 

X = -V^(/f)z = -^z. 

PiK) 

Remark 2: It can also be shown, with an analysis similar to that in Appendix A, that the waves decorrelate over 
frequency offsets larger than e. Thus, one can study the energy density resolved over both time and space i.e., the 
space-time Wigner transform. We do not consider this here, as we limit our study to a single frequency. 


8. Illustration of the scattering effects in statistically isotropic media 

In this section we give numerical illustrations of the net scattering effects in isotropic media. We already stated 
that in transverse isotropic media, where IR depends only on \x\ and z, the matrix Q)*:) defined in (107) is diagonal, 
and depends only on |#t|. We also defined in (28) the scattering mean free paths S(k) and S^iK), which depend only 
on |ic|. The illustrations in this section are for statistically isotropic media, where fR(£) = IRisod^l) and as shown in 
Appendix B, the real part of matrix Q(ic) is a multiple of the identity. Therefore S(k) - S-'-(k), meaning that both 
components of the waves lose coherence on the same l/cl-dependent range scale. 

In Figure 2 we plot the scattering mean free paths for a Gaussian autocorrelation lR(a?) = exp(-|;^p/2), and for 
three different values of y = Aji equal to 27r/50, Injll and 27r/10. We observe that the scattering mean free paths 
decrease with y, meaning that higher frequency waves lose coherence faster. The scattering mean free paths also 
decrease monotonically with |ic|. We may understand this intuitively by noting that at a given range z, the length of 
the path of the waves with transverse vector l/tl is zIP(k). This increases with |*:|, and the waves lose coherence faster 
because they travel a longer distance in the random medium. The monotonicity of the scattering mean free paths can 
also be seen clearly in Section 9.1, where we study the high-frequency limit y —> 0. 

To illustrate the evolution of the coherence matrix ‘P(k,z), satisfying the evolution equation (110), we consider 
first the case of isotropic initial conditions, where Ao(k) depends only on l/cl, and we excite either the transverse 
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Figure 2: Illustration of the scattering mean free paths in an isotropic medium with Gaussian autocon'elation, for three different values of y = Ajt, 
as indicated in the legend. The abscissa is the magnitude |/f| of the transverse wave number. 


magnetic or electric waves. From (63) we see that this corresponds to a current source with - 0 and J along k or 
K^. In the spatial coordinates, the latter means that J is the gradient or the perpendicular gradient of a scalar valued 
function. When the initial conditions are isotropic, the diagonal and off-diagonal terms of the coherence matrix z) 
decouple. Indeed, there is no coupling in the hrst two terms of equation (110), because Q{k) is diagonal. The integral 
in (110) can be analyzed with a Picard iteration. For each iterate we obtain by explicit calculation, using integration in 
polar coordinates, that the coupling terms of the kernel vanish, and that the result preserves the isotropic dependence 
on K. We obtain two autonomous systems of equations for the power vector 


/Tn(*:,z)\ 

\y22(K,z)j 


— lim 


( E[|fl^(^,z)l2] \ 


(115) 


and the off-diagonal terms 

yi2(K,z) - limE[a'^(#i:,z)a'^’-L(#i:,z)]. (116) 

This means in particular that if the wave is transverse electric or magnetic initially, so that 0) = 0, the coherence 
matrix y{K, z) remains diagonal for all z. 

We illustrate the transfer of power between the modes due to scattering by displaying in Figure 3 the z-evolution 
of the power vector (115) for the initial condition 


/Tii(*:,z)\ ^ /|ao(*t)P'l 
\J’22(*t,z)/ \|flo(*t)P/’ 


\ao{K)\^ = exp 



a^(K) = 0 . 


The electromagnetic plane wave at wave vector k is transverse magnetic initially and we have in terms of the Stokes 
vector dehned in ( 1 8) 

e(^,0)-expj-^j(l, 1,0,0). 

In Figure 3 we show the results for y - InjlQ, Jj - y and = 3 . 57 . We plot Tn and CP 22 as functions of l/tl, for 
z = 0, z = 2 X 10“^ and z = 5 x 10“^. Note from Figure 2 that the scattering mean free paths at x = 0 are approximately 
1.3 X 10“^. Thus, in the last line of the plots in Figure 3 the waves have traveled approximately five scattering mean 
free paths. We observe from the hgure that energy from the transverse magnetic waves, which are the only ones 
excited initially, is transmitted to the transverse electric modes. This energy transfer is faster at smaller |*:|, so the 
modes have equal power at hve scattering mean free paths in the case 7 ^ = 7 , whereas for 7 ^ = 3 . 57 , the transverse 
magnetic modes still carry more energy. Recalling the Stokes vector (18), we note that when the transverse electric 
and magnetic waves have equal power, we have &{k, z) = 2 T’i 1 (At, z) (1 , 0 , 0 , 0 ) since CP 12 = 0 in the simulation, and the 
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Figure 3: Display of CP] i {k, z) (full blue line) and 1P22(*^, z) (dashed red line) as a function of |*:|, for three different values of z: Top line is for z = 0, 
middle line is for z = 2 X 10“^ and the bottom line for z = 5 X 10“^. We take y = 2n/50, Jj = y in the left column and jj = S.Sy in the right 
column. 


electromagnetic plane wave at wave vector k is unpolarized. Aside from the transfer of energy between the transverse 
electric and magnetic components of the waves, we note in Figure 3 the diffusion of energy at higher |#t|. For example, 
in the case of - y, the support of the energy is at l/tl <0.1 initially (top left plot) but it extends to l/tl < 0.2 at five 
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Figure 4: Plot of the coefficient of power distribution Cp(z). At z = 0 all the power is in the transverse magnetic mode. As z grows the other mode 
gains power, and eventually the energy becomes equally distributed between the modes. 


scattering mean free paths (bottom left plot). 

Figure 4 is another illustration of the effect of the initial conditions on the power exchange between the modes. 
We consider there, in addition to = y and 3.5y, the case of an even larger y^, equal to 7y. We plot in Figure 4 
the coefficient of power distribution C^(z), calculated as the difference of the mode powers, normalized by the total 
power, which is constant in z. 


!\^\^ydK[yniK,z)-y22{K,z)\ 

- ■ 

{^\^xdK[yn{K,z) + y22{li,z)\ 


(117) 


Consistent with the left bottom plot in Figure 3, is almost zero at z = 5 x 10“^ for y^ = y. However, when jj - ly 
i.e., when the diameter X of the source is 7 times smaller and the opening angle of the emitted wave cone is 7 times 
larger, is still large, so the waves still remember the initial excitation. We terminate the plot for y^ = 7y at z < 0.04, 
where the opening angle of the cone comes close to 180 degrees and the evanescent waves start to get excited. 

The results in Figure 5 display the z-evolution of the entries CP;;(*:, z) of the coherence matrix, for the anisotropic 
initial condition 


\ao{K)f- = exp 


1 

IKX 

1 

(Kl +K2'^ 

2 

1 0.03 / 

2 

1 0.1 / 


a^{K) = 0 . 


Again, we observe the transfer of energy from the transverse magnetic components to the transverse electric ones, 
which is more efficient at l/cl « 0. We also note the diffusion of energy to waves with larger wave vectors, which 
means physically that the focused (beam-like) field emitted by the source spreads out as it propagates along z. The 
off-diagonal components of the coherence matrix are no longer zero, due to the anisotropic initial condition. However, 
the solution approaches an isotropic one as z grows, and the off-diagonal components decay. 


9. High-frequency analysis 

In this section we analyze in detail the net scattering effects of the random medium by considering the high- 
frequency limit y —> 0. We begin with the quantification of the scattering mean free paths and then analyze the 
coherence matrix. 
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Figure 5: Display of Tj i (*:, z) (top line), J’ 22 ('t, z) (middle line) and |lPi 2 (*', z)| (bottom line) as a function of ac = (k[ , K 2 ), for y = 2n/50 and three 
different values of z: Left column z = 0, middle column z = 5 X 10“^ and right column z = 2.5 X 10“^. 


9.1. Scattering mean free paths 

In transverse isotropic random media the matrix Q(/f) is diagonal. We are interested in its real part which dehnes 
the scattering mean free paths (28). We have 


and 


lAfy2 

Re[Qn(.)]=-^ 
Re[Q=WJ = -j^ 


f 


(Ijrf 


Jl 


k{K - k') 




r d(kK')~ 

'k(i(- P)) 

L<1 (2^)' 

i r 




(118) 


(119) 


where we recall that ic- k' - (k- k',I3(k) -I3{k')). Let us change variables k' - k - yu, and expand in powers of y 

(1(k)-P(k - yu) = yu ■ Vf(K) + O(y^) = ^ + O(y^). (120) 
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For the diagonal entries of F we have 


Tji(K-yu,K)- - y — - ¥ 0(y^), i - 1,2, ( 121 ) 

and for the off-diagonal ones 

^ 2 i{k- 7 U,k) = y — -+0{y^), (122) 

\K\\K-yu\ 

and similar for F 12 , with a negative sign. Substituting into (1 18-119) and then into the definition (28) of the scattering 
mean free paths, we obtain using the identity 


that 


,0-(u,u-V/3(k)) 


f du 5i{ku, ku ■ Vj8(/f)) = r du f dr f d^ fR(r, ip}e 

Jb2 Jr2 Jr2 J-oo 

= f dr f d^Ji(r,0(27Tfd[k(r + {VP(K))] 

kJ — OO 


(27ry 

k2 


%j —00 


y 6 (^) 


‘ 5 (^) = Z 








+ 0{y\ 


S^{K)-S{K)^0{y\ 


(123) 


This result shows that when 7 « 1, the scattering mean free paths are almost the same. They are proportional to y and 
decrease as the negative power of 2 with the frequency i.e., the wave number k. That the waves at higher frequency 
lose coherence faster, is expected, as observed in Figure 2. 

Remark 3: If fR(;2) depends only on |^|, i.e. 3?(;2) = !Kiso(|a?|), then we find from (B.21) in Appendix B that 
Re[Q(A:)] is proportional to the 2x2 identity matrix, and 


r 


S{k) (*:) = 77 


ApiK) 




+ O(y^). 


This is in agreement with (123) because 


01 




^ 2 

IPH^) ^ 


— 


(IL 

'\m 


Remark 4: These results hold also when |iit| < y, and we have 

8 


>S(yK) = 




0{y ), S (yK) - Siyic) = 0(y ). 


(124) 


(125) 


This shows that in the high-frequency regime, in the dimensional variables, the scattering mean free paths for small 
wave vectors are equal to S/[a^k^{ fR(0, ^)^f^]. These are the scattering mean free paths obtained in the white- 

noise paraxial regime as described in Appendix C, equation (C. 6 ). The agreement is expected because when y —» 0 
we look at a narrow cone beam propagating through a random medium, as in the paraxial regime. 


9.2. Evolution of the coherence matrix 

In this section we consider the coherence matrix in the high-frequency regime, denoted by 


T’hf(*t,z) = \imy‘^OP{yK,yz), 

y -^0 


(126) 
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and compare it to the coherence matrix in the paraxial regime. We rescaled z by y because we know from (125) 
that the waves lose coherence over such propagation distances, and we see below that energy exchange between the 
transverse electric and magnetic modes occurs at the same scale. We also rescaled *: by y because the initial condition 
generates such wave vectors when y^ is of the same order as y, as we assume here. Taking the limit of (110) as y ^ 0, 
and using 


limr(y#t, y/t') = rhf(/!:,/f') 

y^O 


1 Ik k' 


(127) 


and 


limyQ(y*:) = Qhf(*t) 

y -^0 




f 

Jr 


du 


8 Jb 2 (2;r)^ 


Jliku, 0 ) 


(128) 


we find that T^hf (*^, z) satisfies 


d^TMiK,z) =2Qm (#t) J’hf (*:, z) + 



dikiJ) 


rhf(A:, /c')J’hf('t', z)rhf(#t', K)Ji{k(K - iJ), 0), 


(129) 


starting from J’m(*’, 0 ) = ./lhf,o(*^)Ahf,o(#t)^ with 


'^hf,o(^) 


/ ^hf,o(^)\ 


flhf,o(^) — T| I ' 

2\k\ 





(130) 


and Kj — limy_>o7//7 = Oil). More explicitly, equation (129) written componentwise is 


5 zThf,ll(>f, 2 ) = ^ j 
(a: ■ k')- 

^ (Z, \1 

B 2 ( 22 r )2 

^yM.niK',z) - {k ■ k')(k^ 

[ - J’hf,ii(A:,z)+ 

■ *t')(J’hf,i 2 ('i:',z) + J’hf. 2 i('t',z)) + ■ *:')^lPhf, 22 ('t',z) 

{k ■ K'y 

^ (Z, ',2 

b 2 ( 22 r )2 

(a:-^ ■ K’){ 7 M^n(K',z) - 7 \ 

l/tpl/f'p 

[ - 7 M.niK.,z)+ 

il, 2 liK',z)) + {k ■ #t')^lPhf, 12 (*t',z) - ■ *: 0 ^ 1 Phf, 2 l('t',z) 

5 zThf, 2 l(*t,z) = ^ J 

(k ■ K'y 

r 2 

B 2 ( 22 r )2 

■ K’){ 7 M,l\iK',z) - 7 \ 

l/tpl/t'p 

[ - J’hf, 2 l(A:,z)+ 

if, 22 (*t',z)) - (^-^ ■ K'fyu.uili^z) + {k ■ K')^ 7 M, 2 \iK',z) 

5 zJ’hf, 22 (*t,z) = ^ J 

(k-^ ■ k' 

^ O ',2 

B 2 ( 22 r )2 

)^ 5 ’hf.ii('f',z) + ('f ■*:')(*■ 

l/tpl/t'P 

[ - J’hf, 22 (*^, 2 )+ 

■ *:')( 5 ’hf,i 2 (*:',z) + 5 ’hf. 2 i('f',z)) + {k ■ *:')^yhf, 22 (*’',z) 




This system can be simplified by changing of coordinates, as follows. Let us introduce (for k - (a:i, /C 2 )) 


Kxa^iK,z)-K2a'^’^iK,z) 
ay{K,z) - -—- 


alia, z) 


K2a^iK,z) + Kxa’^’-^iK,z) 


(131) 


and note from (94) that they are the Fourier coefficients of the electric field on the xi and X 2 axes. Define also the 
coherence matrix 


T’hf('t,z) = limy'‘lP(y#i:,yz), 

y^O 


lPiK,z) = limE 


/ai(*:, 2 )\ 


/a^(#t,z)\^ 

\aliK,z)j 


(132) 
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which is related to CPhf (*■, z) as 


~ , , K\7i,{^n{K,z) +Kl7h{,22(‘(,z)-KxK2(y^f^2(K,z) + 7hf,2\{X-^^)) 

7^f,n{K,z) = -!-^- 

kP 

~ KiK2{7M,l\{K;Z)-7bf,22{li,z)) +K\7hf,n{K,z)-Kl7),i2\{K,z) 

yhf,i2(^,z) =-^- 

kP 

yM,2i(^,z) =-nr^-^- 

kP 

yhf,22(^,z) = ^^- 


Elementary calculations then show that iPhf (^, z) satisfies 

_ 1^ (y^ r* cJ(Jck'^ r _ _ 1 

d^7hf(K,z) = —— I 7{k(K - *:'),0)f - J’hf(*:,z) + J’hf(#t',z)], (133) 

4 Jk 2 (27r)^ L j 

starting from 7hf(K, 0) = ^hf,o(*^)Ahf,o(#t)^ with 





(134) 


This equation agrees to the one found in the white noise paraxial regime, as described in Appendix C, equation (C.7). 
It says that in this high frequency scaling there is no polarization exchange. If the electric field is along the xi axis 
initially, it remains along the xi axis for distances that are comparable to the scattering mean free path. 

Remark 5: The polarization conservation described here happens only in the high-frequency regime, where the 
transport equation reduces to (133). It does not hold in general, and equation (110) exhibits polarization exchange. 


10. Summary 

We presented a detailed analysis of cumulative scattering effects on electromagnetic waves that propagate long 
distances in random media with small fluctuations of the wave speed. The results are derived from first principles, by 
studying Maxwell’s equations with random electric permittivity, in a regime of separation of scales modeled by two 
parameters. The first is e, the ratio of the wavelength and the distance of propagation, and the second is y, the ratio of 
the wavelength and correlation length of the random fluctuations, which is similar to the spatial support of the source. 
By controlling y, we ensure that the source emits a cone wave beam which propagates along a preferred direction, 
called range. Depending on y, the opening angle of the cone may be much larger than in paraxial regimes. However, 
the angle is smaller than 180 degrees, so that evanescent and backscattered waves can be neglected. 

The analysis of the solution of Maxwell’s equation is in the long range limit e —> 0. It involves the decomposition 
of the waves in transverse electric and magnetic plane wave components, whose amplitudes are random fields. They 
satisfy a stochastic system of differential equations driven by the random fluctuations of the electric permittivity. 
These equations model the evolution in range of the random amplitudes, starting from the initial values determined by 
the source of excitation. The e —» 0 limit of the solution of the system of stochastic differential equations is obtained 
with the Markov limit theorem. 

We study in detail the first two statistical moments of the limit amplitudes. Their expectation decays exponentially 
with range, on length scales called scattering mean free paths. This is the manifestation of the randomization of the 
waves, due to scattering in the random medium. The second moments of the amplitudes describe the decorrelation 
of the waves over directions, and the transport of energy. Of particular interest is the energy density (mean Wigner 
transform), which characterizes the state of polarization of the waves and the diffusion of energy over directions. The 
transport equations with polarization satisfied by the energy density are derived with the Markov limit theorem, and 
the result is related to the radiative transport theory. In the high frequency limit y —> 0 the equations simplify, and 
are related to those satisfied by the Wigner transform of the solution of the paraxial wave equation in random media. 
We quantify the transfer of energy between the components of the wave field, and illustrate the results with numerical 
simulations. 
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Appendix A. The Markov limit 

Let 0 be an open set in and G{0, K^) the space of infinitely differentiable functions with compact support. We 
consider the process in C([0, L], D'), the solution of 

(A.1) 

uz yje \e e) \e e) 

where ^') and are random linear operators from £)' to £)'. We assume that the mappings f ^ 

and ^ ^ (') are stationary and possess strong ergodic properties, and that FiF 0) has mean zero. Moreover, the 

mappings 0 F(F C) ^nd 0 —> @(F 0) are periodic. 

We are interested in particular in equation (102), where the process is in £)'(0, M^), defined by 

' Re(fl^(*:, z))' 

, for a: e 0 = {x € l/cl < 1). 




Im(fl^(*:, z)) 
Re(a"’-"(x,z) 
Im(a'^’-^(A:,z)/ 


(A.2) 


The operator Fi0 C) is 

{F(Fnx,<i>) = y f d^[F(Fn^]jmj('^) ^ f dKm- f 

Jo Jo Jo 


), 


(A.3) 


for ^ e £)(0,K'*) with components and X 6 £)'(0,IR"^) with components A/. The kernel matrix in 

(A.3) is given by 


F = 


^ 11 

-F 

11 

^ 12 

-T \ 
^ 12 

^ 11 

“^11 

rri 
^ 12 

^ 12 

^21 

-F 

“^21 

^22 

"^22 


F 

rri 

^22 

rrr 

^22 ^ 


in terms of 




L 2(27ryy^ 


^ , r 


7 

k{K - k') 

7 


rCjFJiK,K’,nl 

y^F^;{K,K',a\ 


(A.4) 

(A.5) 

(A.6) 


L 2{2TJp-y^ 

where we recall from (76)-(80) the expression of F‘‘f{K, k' ,0). The adjoint operator F*{0 C) is defined by 

{F(FOX,<l>)^{X,F\C,0<P) 

for 0 6 £)(0, K"^) and X e iD'(0, M'^), and has kernel F*(a:, k',0 0) - F'^{k',k, 0 0). Similar expressions hold for g. 

To obtain the Markov limit we use the results in [12]. The interested reader can find a self-contained introduction 
to such limit theorems in [20, Chapter 6] or in [21, Appendix A]. We get that X^(z) converges weakly in C([0, L], 2)') 
to X(z), which is the solution of a martingale problem with generator £. defined by: 

2:/«X,0))= r dC\im]-{ dhF[{X,F\O,h)<l>){X,F\0^ + h)<p)]f"{{X,<p)) 

Jo 4 Jo 


f Jfhmi f dhE[(X,F*(O,h)F'(0^ + h),p)]f((X,,p)) 
Jo 4 Jo 

1 

- dhE[(X,0*(O,hm]f'((X,<fi}), 

^ Jo 


+ lim — 

Z—»oo Z 


(A.7) 
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for any X € £)'(0,K"^), 0 e £)(0,K"^), and smooth / ; K —» M. This means that, for any 0 6 £)(0,IR"^) and smooth 
function f :M. M., the real-valued process 


f({X(z), 0» - r dz' m{X{z), (P)) 

Jo 

is a martingale. More generally, if n 6 N, 0i,..., 6 D(0, K"^), and / : K" ^ K is a smooth function, then 

f((X(z), (X(z), 0„» - r dz' £'fi{X{z'), 0i),..., <X(z'), (/>„)) 

Jo 


is a martingale, where 

rf({x,,p,{x,ci>„)) 

Jv r 


1 I 

hm - J dhE[{X,r*(Q,{X ,^ h- /!)0,> ] d]i 
1 

d( lim - dhE[ lx, r*(0, h)r\(, ^ + h)^,) ] 5/ 

z^oo Z Jo ^ ■ 

J Jim i dhE[ {X, r (0, h)<Pj) ]dj MX, 4,,{X, (P ,)). 


n 

5 J. 


To calculate the hrst moment of the limit process X(z), let n = 1 and f(y) = y in (A.8). We find that 

dE[{X{z),(P)] 


dz 


= e[ {x(z), T(p) ] + e[ {x(z), g(p)]. 


where 


r = 


= f d4 lira - f dhE['F*(0,h)'FM{ + 
Jo z Jo 

f dhE[g*(0,h)]. 

Jo 


1 

g = lim — 
^ Z^oc Z 


This shows that 


X(z) = E[X(z)], 

satisfies a closed system of ordinary differential equations 


d{X{z),<p) ^ ^ {X{z),g(p), 


dz 


or, equivalently in £)', 


'^-^^TX(z) + gX(z), 

dz 


where T, resp. g, is the adjoint of T , resp. g . 

Recalling from (A.3-A.6) the expression of the kernel F^(/c', k, g g) of T^ig C), we obtain 

1 

— I c/h E[F;^/!)F^j(/!:",*:, 0 ,/!)], 




(A.8) 


(A.9) 
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for j,l - 1,..., 4. For instance. 


Tix{ic,k')- f cIk” f d(\\m— f c/ZrE[9^j^ 0,/;)] 

Jo Jo Z Jo 

- f dtc” f d(\\m— f dhE[j\^(/(',K",^,^ + h)j\i(/(",/c,0,h)] 

Jo Jo Z Jo 

+ r dtc” f cZ^ lim ^ r 

Jo Jo Z Jo 

- f dK” f £/^ lim ^ r £//!E[Ji 2 (#t',/!:",^,^ 0,Zi)], 

Jo Jo ' Z Jo 


and using (A.5-A.6), we get 


r 


( Jc^ C* 1 

sw) v""!, 


k(/c' - k") 
7 


jnv 


k(K" - k) 
7 


,0 


[F‘“‘{k', k",^ + h)F‘^‘'{K'',K, Zl)] 1, 


Moreover, using the identity 

{k{K'-K") 


E 


7 


^7^ y 


k(K" - k) 
7 


,0 


= 1^1 


k{K - k") 




derived from the dehnition of the autocorrelation with straightforward algebraic manipulations, and obtaining from 
(76) that 




we get 


x[r'’“(A:, k")T‘'%k'', *:)] 1 j 6{k - a:')}. 

The expressions of the other components of 'Fji(k,k’) are of the same type. Similarly, we calculate Q , and 
substituting into (A.9) we obtain the explicit expression of the differential equations satished by the mean wave 
amplitudes. This is equation (106), written in complex form. 

The calculation of the second moments is similar, by letting « = 1 and f(y) = in (A. 8), and carrying the lengthy 
calculations. 


Appendix B. Connection to the radiative transport theory 

To connect our transport equations (112) (see also (110)) to the radiative transport theory in [13, 9, 10], we adhere 
to the notation in [9]. First note that in [9] the random medium is assumed to be statistically isotropic, i.e. !li(®) 
depends only on |a?|. 

For any 3C e consider the orthonormal basis {;?*°'(9C), of dehned by 

^^^)=^ = ^(3C,9C,), /2)(3^) = (^,o), (B.l) 

|9C| |3C| |9C|\ W / \|3C| / 

and satisfying 

^0)x^l)^^2)^ ^0) X ^2) ^ _^1)^ (2 2) 
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In [9] the vector is denoted by k. We call it to avoid confusion with the wavenumber k. Note that when 

DC - kit, with = 1, then = k, and the vectors and are the same as it and defined in (17). 

Following the notation in [9], we define 

//lie, iS) = X ], j = 1,2, (B.3) 

where we use a different constant of proportionality than in [9], to simplify the relation (B.5). The Wigner transform 
W(IX, £) is the 2x2 matrix with components 


W,„(3C, 


l,q^l,2. 


(B.4) 


where the bar denotes complex conjugate. Substituting the wave decompositions (94-95) into (B.3), and keeping only 
the forward propagating modes, we obtain after some algebraic manipulations that 


W(3C/k, x, z), 

with W the Wigner transform in (1 1 1), satisfying (112). 

The transport equation in [9] is 

for the dispersion relation 

co(DC) = cJ3C|. 

The integral kernel in the right hand side of (B.6) is the differential scattering cross section, defined by [9] 


^ TTc^k^a^ ~ I DC - DC 

cr(3C,fK: )[W(3C ,®)] = . 3^ 


2(27r)^y^ 


with the 2x2 matrix 


T(^,Dt') ^ (Y\Dt) ■ 

The scalar S(1K) is the total scattering cross section, which is shown in [9] to satisfy 

S(1K)I = J dXcr(DC,DC')[I]. 


(B.5) 


(B.6) 


5[aj(DC) - lo{DC )]T(3C, DC )W(fK;, x)T{DC , DC), (B.7) 


(B.8) 


Since (B.5) gives that W(9C, a?) is supported at vectors DC of the form DC = kit, with it = {k,I3{k)), we see that the 
operator on the left hand side of (B.6) is given by 


Vw(DC) ■ = Co— ■ Vai = Co[k ■ Va, -HP(/t)5z] = CoP(/f)[Pz - Vp(*:) ■ Va,]. (B.9) 

|dC| 

Again, using (B.5), we see that the integral kernel in the right hand side of (B.6) is supported at vectors DC = kit, 
with it = (k',/3{k')), so the Dirac distribution in (B.7) is 

d[cu(DC) - caikit)] = d[co ^IDCp + |DC,p - c^k] = ■ (B-IO) 

Thus, the integral in (B.6) is supported at DC = kit, with it = (k,I3(k)). For such vectors, the matrix T equals 

Tikit, kit) = yll3{K)li(K')T{K, k’). (B. 11) 
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From (B.5) and (B.9), we obtain that 


Vw(3C) ■ V£W(9C, £) = Cod[X, - k^(DC/k)][d, - VpiXilk) ■ V^]W(aC/A:, x,z\ (B.12) 


from (B.5), (B.7), and (B.lO-B.ll), we obtain that 


/ 




r k'^rp- r 

V Jk'isi 


d(kK') - i Xi-kK' kiJiCKIk) -P(k')) 

(27r)2 \ r ’ r 


X 


T(JC/k, k')W{k', X, z)r(K\ DC/k). 


(B.13) 


We also find from (B.8) that 


S(aC)I = 


clkV 

4(27r)2y^ 

or equivalently, for DC = k{K,f5(K)), 

Cok^P{K)a^ 


J djt'6[a>(jt') - 


OC-X 

7 


T(3C,3C)T(3C ,3C), 


so that 


2(3C)I = 


s(ac)W(ac, x) = 


V Jk'i 

Cgk^a^ 

4y3 


r d{kK' 

Lki 


)~(k{K-K’) k(j3{K)-f5(K’)) 


r(*:, 


6[X,-kf5(X/k)] f 

Jk'l<i 


d(k/c')-/X-k/d kipiXlk) -I3 (k')) 

7 ^ 


(B.14) 


(B.15) 


(B.16) 


TiXlk, k')Y{k’, Xlk)W(Xlk, x,z). 

Finally, using the transport equation (112) satisfied by W)*:, x,z), and the relation (B.5), we obtain that 

iaj(X) ■ ViiW(lK, JdX o-(X, ljc')W(ljc', x) + Coli(Xlk)[Q{Xlk)W{X, x) + W(3C, x)Q{Xlk)^\ (B.17) 


This is similar to the transport equation (B.6), except that we do not have the scalar valued total scattering cross 
section S(3C) multiplying W(3C, :2), but a linear operator -CoyS(3C/A:)|^Q(3C/A:)W(3C, ;^) + W(3C, £)Q(3C/A;)’j acting 
on W(IX, a?). The results would be exactly the same if Q)*:) were a multiple of the identity of the form 


-2c„p{KmK) = S(DC)I 


for X - k(K,/3(K)). This turns out to be the case in the high-frequency limit y ^ 0, as explained in Section 8. It is 
also the case for the real part of Q(/i:) in general (provided that is isotropic, as assumed throughout this appendix), 
as we now show. 

For 6 such that l/tl < 1, we have by (108) 


Re(Q(^)) = - 


8 r^ J\K'\ 


f d(kid)~ k(K- k) k(Ji{K) - P{k))\ , 

77 - 5 - 3 ? -,- r(/i:,ii:)r(/t,#t). 

(2;r)2 \ y 7 


By comparing with (B.15) and denoting X - k(K,/3(K)), this shows that 

-2c„A^)Re(Q(^)) = 5:(lj^)I, 


(B.18) 


(B.19) 


which is a multiple of the identity matrix. This also follows by inspection from (B.18), rewritten as 

kW 


Re(Q(^)) = -- 


H2^)^7W M'\=k 


JlJC l=k 


dS(X )X 


'’x-ot^ 


7 


T(X,X)T(X ,X). 


(B.20) 
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Indeed, if we parameterize 


then we have 


where 


I sin 9 cos (j) 
ysin6>sin0 


'sin 0 cos 0' 

II 

'sin 9' cos (p'^ 

sin0sin<;ii 

sin 9' sin p' 

cos 9 


cos 9' 




l-« 


1 - k'; 


\DC-DC\^ = 2k^il-K'^), 



(k-"\ 

'‘i 


sin 9' sin(0' - (p) 

K = 

k” 


cos 0 sin 0' cos(0' - (p) + sin 0 cos 0' 


k" 


^sin 0 sin 0' cos(0' — (p) + cos 0 cos 0' ^ 


If we introduce the rotation matrix 


U. 


- sin 

cos 0 cos (ft 
sin 6* cos 0 


cos (f) 
cos 0 sin 0 
sin 9 sin ((> 


0 ' 
- sin0 
cos 9 ^ 


then 

U^IK' = kK. 

By carrying out the change of variable k' - 1]^% jk in (B.20), we find that since [R(i?) = !Kiso(|^|), 


Re(Q(^)) = 




UlnYy^PiK) 


f 


dS(k")0lis 


V2(r 


•^r") 


^2 




-<4 

1 - Kf 


(B.21) 


By changing into and carrying the integration over the unit sphere, we obtain that Re(Q(iif)) is 

proportional to the identity matrix. 


Appendix C. Connection to the white-noise paraxial theory 

The white-noise paraxial theory is valid when the four length scales satisfy 

A-^e-X-^L, (C.l) 

with AL ~ X^, which corresponds to a Fresnel number of order one. The separation of scales can then be characterized 
by a unique dimensionless parameter e « 1 such that 


e - 


A 

V 



(C.2) 


We again let L be the reference length scale and introduce the scaled length variables 

x'^xKeil), z'^zll, L'^LIL, X'^XHeiL). (C.3) 


The scaled wavenumber is k' - kLe - 2n. We finally assume that the standard deviation a of the fluctuations of 
the random medium in (5) is small, of order and we let a' - be the normalized standard deviation of 

the fluctuations. If we drop the primes and consider the asymptotic regime e —» 0, then the electric field W has the 
form [8] 




/Ep^(X, Z)\ 

I 0 / 




r Tp^Yx,x',z)Eo{x')dx', 

Jr2 


(C.4) 
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where the scalar transmission kernel satisfies the Ito-Schrodinger equation 


I IK 

dTp^(x,x',z) = —A^,Tp^{x,x',z)dz+ —T^^{x,x',z)odB{x',z), (C.5) 

starting from T-p^^ix, x',z -0) - 6{x - x'). Here B(x,z) is a Brownian field with mean zero and covariance function 

¥.[B{x,z)B{x',z')\ = rm.n{z,z')a^t J" 

and o stands for the Stratonovich integral. Using Ito’s formula, the coherent field Ap.^-(x,z) - E[£^pai(£c,z)] satis¬ 
fies [8] 


^z-^par(^;^) — o -^par(^5^)5 

2k .Spar 

which shows that it decays exponentially with z on the scale 

g 

~ ^2^2 J_” :R(0,Od^' 

The decay length .Spar corresponds to the scattering mean free path S(k) defined by (123) in the high-frequency regime 
y —» 0 with |*:| = 0(y). 

The Wigner transform is 

Wpar(*:,at,z) = J^dye'‘"‘">^E^Ep^(x - ^,z)Ep^,(x+ 

" £2 £$ +\\^ ^par - I)’ ^) ]■ 

Using Ito’s formula it is shown in [8] to satisfy the transport equation 

5^Wpar + ■ Va^Wpar = f - k ’), 0)[Wpar(^') - Wpar(^)], (C.7) 

4 Jb 2 \2ny 

which corresponds to (133) in the high-frequency regime y —> 0. 
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